#########################################################################   
#                           INFO                                        #   
#########################################################################  

# PROJECT: Gender Stereotypes, Political Leadership, and Voting Behavior in Tunisia
# PURPOSE: Correlates of patriarchal attitudes & residuals robustness check
# CREATED: September 2019 by Alexandra Blackman & Marlette Jackson
# INPUTS:  tun_bjka.xlsx; tun_yg.xlsx
# OUTPUTS: Correlates of Patriarchal Index table; robustness check with econ priority 
#           respondents; robustness check with residuals

#########################################################################   
#                         SETUP                                         #   
######################################################################### 

rm(list = ls()) 
setwd("~/Dropbox/Voters and Female Politicians/Political Behavior")


#########################################################################   
#                        LOAD PACKAGES                                  #   
######################################################################### 

need <- c("doBy","foreign", "lubridate","extrafont", "car", "plyr",
          "reshape", "openxlsx", "tidyr", "dplyr", "ggplot2", "ggthemes",
          "RColorBrewer", "rms", "lmtest", "sandwich", "msm",
          "stargazer","texreg","cowplot","plm","margins","prediction",
          "gplots","ggpubr", "sjmisc","sjPlot","gtable","cjoint","cregg",
          "scales","multiwayvcov") # packages needed
have <- need %in% rownames(installed.packages()) 
if(any(!have)) install.packages(need[!have]) 
invisible(lapply(need, library, character.only=T)) 

#########################################################################   
#                 LOAD CLUSTERING, F-TEST FUNCTIONS                     #   
#########################################################################   

super.cluster.fun <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}

f.test <- function(model, cluster){
  vcovCL<-cluster.vcov(model, cluster)
  ftest <- waldtest(model, vcov = vcovCL, test = "F")
  return(ftest)
}

#########################################################################   
#                       LOAD DATA                                       #   
#########################################################################   

## LOAD BJKA HOUSEHOLD SURVEY (JULY 2017)
tun_bjka <- read.xlsx("tun_bjka.xlsx")

## LOAD YOUGOV ONLINE SURVEY (APRIL 2019)
tun_yg <- read.xlsx("tun_yg.xlsx")


#########################################################################   
#         ANALYSIS OF CORRELATES OF PATRIARCHAL ATTITUDES               #   
#########################################################################

# Q: What are the individual level correlates of patriarchal attitudes?
table(tun_bjka$patr_index)
table(tun_yg$patr_index)

# Basic model w/ delegation FE
patr_bjka1 <- lm(patr_index ~ resp_fem + age + del_en, data = tun_bjka)
summary(patr_bjka1)

# Basic model w/ city FE
patr_yg1 <- lm(patr_index ~ resp_fem + age  + city_en, data = tun_yg)
summary(patr_yg1)

# Controlling for education, marital status, voting, employment
patr_bjka2 <- lm(patr_index ~ resp_fem + age + educ_high +
                   resp_mar_wid + vote_yes + employ_status + del_en, data = tun_bjka)
summary(patr_bjka2)

# Controlling for education, marital status, voting, employment
patr_yg2 <- lm(patr_index ~ resp_fem + age + educ_high +
                 resp_mar_wid + vote_yes + employ_status + city_en, data = tun_yg)
summary(patr_yg2)

# Controlling for religiosity 
patr_bjka3 <- lm(patr_index ~ resp_fem + age + educ_high +
                   resp_mar_wid + vote_yes + employ_status + 
                   relig_great_imp + attend_daily + del_en, data = tun_bjka)
summary(patr_bjka3)

# Controlling for religiosity 
patr_yg3 <- lm(patr_index ~ resp_fem + age + educ_high +
                 resp_mar_wid + vote_yes + employ_status + 
                 relig_great_imp + attend_daily + city_en, data = tun_yg)
summary(patr_yg3)

## Robust standard errors
patr_bjka2_ro <- coeftest(patr_bjka2, vcov = vcovHC(patr_bjka2, type="HC3"))
patr_bjka3_ro <- coeftest(patr_bjka3, vcov = vcovHC(patr_bjka3, type="HC3"))
patr_yg2_ro <- coeftest(patr_yg2, vcov = vcovHC(patr_yg2, type="HC3"))
patr_yg3_ro <- coeftest(patr_yg3, vcov = vcovHC(patr_yg3, type="HC3"))

# Output table format
stargazer(patr_bjka2, patr_bjka3,patr_yg2, patr_yg3,
          title = "Correlates of Patriarchal Attitudes",
          dep.var.labels   = c("Patriarchal Attitudes Index"),
          covariate.labels = c("Female", "Age: 30 to 39","Age: 40+", "Secondary Education+", 
                               "Married/Widowed", "Voted in most recent election", 
                               "Unemployed or out of labor force",
                               "Religion is of great importance","Attend prayer daily"),
          omit = c("del_en","city_en"),
          font.size = "small",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE,
          notes = "Robust SEs",
          add.lines = list(c("Delegation/City FE" ,"Yes","Yes","Yes","Yes"),
                           c("Sample" ,"BJKA","BJKA","YouGov","YouGov")))


# Output table format
stargazer(patr_bjka2_ro, patr_bjka3_ro,patr_yg2_ro, patr_yg3_ro, 
          title = "Correlates of Patriarchal Attitudes",
          dep.var.labels   = c("Patriarchal Attitudes Index"),
          covariate.labels = c("Female", "Age: 30 to 39","Age: 40+", "Secondary Education+", 
                               "Married/Widowed", "Voted in most recent election", 
                               "Unemployed or out of labor force",
                               "Religion is of great importance","Attend prayer daily"),
          omit = c("del_en","city_en"),
          font.size = "small",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE,
          notes = "Robust SEs",
          add.lines = list(c("Delegation/City FE" ,"Yes","Yes","Yes","Yes"),
                           c("Sample" ,"BJKA","BJKA","YouGov","YouGov")))


#########################################################################   
#       CORRELATES OF ALWAYS RESPONDERS (BJKA SURVEY ONLY)              #   
#########################################################################

tun_bjka$never_resp <- ifelse(is.na(tun_bjka$choice_1)==T &
                                is.na(tun_bjka$choice_3)==T &
                                is.na(tun_bjka$choice_5)==T &
                                is.na(tun_bjka$choice_7)==T, 1,0)


tun_bjka$sometimes_resp <- ifelse(is.na(tun_bjka$choice_1)==T |
                                    is.na(tun_bjka$choice_3)==T |
                                    is.na(tun_bjka$choice_5)==T |
                                    is.na(tun_bjka$choice_7)==T, 1,0)

tun_bjka$cjoint_response <- ifelse(tun_bjka$sometimes_resp==1,"Sometimes","Always")
tun_bjka$cjoint_response <- ifelse(tun_bjka$never_resp==1,"Never",tun_bjka$cjoint_response)
table(tun_bjka$cjoint_response)

tun_bjka$always_resp <- ifelse(tun_bjka$cjoint_response=="Always",1,0)

# Controlling for education, marital status, voting, employment
always_bjka2 <- lm(always_resp ~ resp_fem + age + educ_high +
                     resp_mar_wid + vote_yes + employ_status + 
                     del_en, data = tun_bjka)
summary(always_bjka2)

# Controlling for religiosity 
always_bjka3 <- lm(always_resp ~ resp_fem + age + educ_high +
                     resp_mar_wid + vote_yes + employ_status + 
                     relig_great_imp + attend_daily + del_en, data = tun_bjka)
summary(always_bjka3)

#########################################################################   
#     ROBUSTNESS CHECK 1: ECON PRIORITY RESPONDENTS ONLY (BJKA only)    #   
#########################################################################

tun_bjka$prior_econ <- ifelse(tun_bjka$priority== "Fight corruption" |
                                tun_bjka$priority== "Improve the economic situation" |
                                tun_bjka$priority== "Job creation", 1,0)


# Subset to just econ priority only
tun_bjka_econ <- tun_bjka[tun_bjka$prior_econ==1,]

#########################################################################   
#       robustess check 1: BJKA vignette analysis                       #   
#########################################################################

m_all <- lm(support ~ tassign_women, data=tun_bjka_econ)
summary(m_all)

m_patr <- lm(support ~ tassign_women, data=tun_bjka_econ[tun_bjka_econ$patriarchal==1,])
summary(m_patr)

m_egal <- lm(support ~ tassign_women, data=tun_bjka_econ[tun_bjka_econ$patriarchal==0,])
summary(m_egal)


# Extract all coefficients 
bjkaFrame <- data.frame(var = rownames(summary(m_all)$coef),
                        coef = summary(m_all)$coef[, 1],
                        se = summary(m_all)$coef[, 2],
                        modelName = "Overall")

patrFrame <- data.frame(var = rownames(summary(m_patr)$coef),
                        coef = summary(m_patr)$coef[, 1],
                        se = summary(m_patr)$coef[, 2],
                        modelName = "Patriarchal")

egalFrame <- data.frame(var = rownames(summary(m_egal)$coef),
                        coef = summary(m_egal)$coef[, 1],
                        se = summary(m_egal)$coef[, 2],
                        modelName = "Egalitarian")

# merge data
allModelFrame <- data.frame(rbind(bjkaFrame, patrFrame, egalFrame))

# drop Intercept estimate
allModelFrame$drop <- ifelse(allModelFrame$var %in% 
                               grep("(Intercept)", allModelFrame$var, value=TRUE), 1, 
                             ifelse(is.na(allModelFrame$var), NA, 0))
allModelFrame <- allModelFrame[allModelFrame$drop==0,]

allModelFrame$var2 <-NA
allModelFrame$var2[allModelFrame$var=="tassign_womensecurity"] <- "Miriam Security"
allModelFrame$var2[allModelFrame$var=="tassign_womentraditional"] <- "Miriam Traditional"
allModelFrame$var2[allModelFrame$var=="tassign_womenwomen"] <- "Miriam Women"


allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Egalitarian",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Patriarchal",
                                               drop=0, var2="Miriam Control")

allModelFrame[nrow(allModelFrame) + 1,] = list(var=NA, coef=0, se=0, modelName="Overall",
                                               drop=0, var2="Miriam Control")

allModelFrame$var2 <- factor(allModelFrame$var2, levels = c("Miriam Control",
                                                            "Miriam Women",
                                                            "Miriam Traditional",
                                                            "Miriam Security"))

allModelFrame$modelName <- factor(allModelFrame$modelName, levels = c("Egalitarian",
                                                                      "Patriarchal",
                                                                      "Overall"))

# Specify the width of confidence intervals
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
vignette_bjka_econ <- ggplot(allModelFrame, aes(x=var2, y=coef, colour=modelName)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  #coord_cartesian(ylim=c(-.5,1)) +
  ylim(-.5,1.2) +
  coord_flip() +
  geom_point(aes(x = var2, y = coef, shape = modelName, colour=modelName),
             size = 3, lwd = 1/2, position = position_dodge(width = 1/2)) +
  geom_errorbar(aes(ymin = coef - se*interval2, ymax = coef + se*interval2), width=.15,
                lwd = .75, position = position_dodge(width = 1/2)) +
  scale_shape_manual(values=c(16, 15, 17),name="Type", guide = guide_legend(reverse = TRUE)) +
  scale_color_manual(values=c('#56B4E9','#E69F00', '#999999'),name="Type", guide = guide_legend(reverse = TRUE)) +
  theme_bw() +
  theme(axis.text.x=element_text(size=12),
        axis.title.x=element_text(size=12),
        axis.text.y=element_text(size=10),
        axis.title.y=element_text(size=12),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank()) +
  xlab("Vignette Treatment") + 
  ylab("Coefficient")

vignette_bjka_econ

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/vignette_bjka_econ.png",
    width = 6, height = 6, units = 'in', res = 300)
vignette_bjka_econ
dev.off()


#########################################################################   
#       robustess check 1: BJKA conjoint analysis                       #   
#########################################################################

tun_bjka_econ$never_resp <- ifelse(is.na(tun_bjka_econ$choice_1)==T &
                                is.na(tun_bjka_econ$choice_3)==T &
                                is.na(tun_bjka_econ$choice_5)==T &
                                is.na(tun_bjka_econ$choice_7)==T, 1,0)

tun_bjka_econ$sometimes_resp <- ifelse(is.na(tun_bjka_econ$choice_1)==T |
                                    is.na(tun_bjka_econ$choice_3)==T |
                                    is.na(tun_bjka_econ$choice_5)==T |
                                    is.na(tun_bjka_econ$choice_7)==T, 1,0)

tun_bjka_econ$cjoint_response <- ifelse(tun_bjka_econ$sometimes_resp==1,"Sometimes","Always")
tun_bjka_econ$cjoint_response <- ifelse(tun_bjka_econ$never_resp==1,"Never",tun_bjka_econ$cjoint_response)
table(tun_bjka_econ$cjoint_response)


tun_bjka_econ_long <- reshape(data = tun_bjka_econ, 
                         idvar = "responseid", 
                         varying = list(Sex=c("sex_1","sex_2","sex_3","sex_4","sex_5","sex_6","sex_7","sex_8"),
                                        Party=c("party_1","party_2","party_3","party_4","party_5","party_6","party_7","party_8"),
                                        Job=c("job_1","job_2","job_3","job_4","job_5","job_6","job_7","job_8"),
                                        Town=c("town_1","town_2","town_3","town_4","town_5","town_6","town_7","town_8"),
                                        Education=c("educ_1","educ_2","educ_3","educ_4","educ_5","educ_6","educ_7","educ_8"),
                                        Age=c("age_1","age_2","age_3","age_4","age_5","age_6","age_7","age_8"),
                                        Policy=c("policy_1","policy_2","policy_3","policy_4","policy_5","policy_6","policy_7","policy_8"),
                                        Experience=c("exper_1","exper_2","exper_3","exper_4","exper_5","exper_6","exper_7","exper_8"),
                                        choice=c("choice_1","choice_2","choice_3","choice_4","choice_5","choice_6","choice_7","choice_8"),
                                        contest=c("contest_1","contest_2","contest_3","contest_4","contest_5","contest_6","contest_7","contest_8")), 
                         direction="long", 
                         v.names = c("Sex","Party","Job","Town","Education","Age","Policy","Experience","choice","contest"),
                         sep="_")

tun_bjka_econ_long_na2 <- tun_bjka_econ_long %>% filter(!is.na(choice)) %>% filter(!is.na(patriarchal)) %>%
  mutate(Sex = factor(Sex, levels=c("Male","Female"))) %>%
  mutate(Party = factor(Party, levels=c("Initiative Party","Afek Tounes","UPL","Democratic Current",
                                        "Current of Love","al-Irada Movement","Popular Front",
                                        "Machrouu Tounes","Ennahda","Nidaa Tounes"))) %>%
  mutate(Job = factor(Job, levels = c("Director", "Lawyer","Teacher","Employee in the private sector","Political Activist","Union Activist"))) %>%
  mutate(Town = factor(Town, levels = c("within 10 km of here","within 20 km of here","within 30 km of here","within 40 km of here"))) %>%
  mutate(Education = factor(Education, levels = c("Bachelors","Masters","Studied outside of Tunisia"))) %>%
  mutate(Age = factor(Age, levels=c("70","60","50","40","30"))) %>%
  mutate(Policy = factor(Policy, levels=c("Improve security","Improve women's rights","Fight corruption","Improve the economic situation"))) %>%
  mutate(Experience = factor(Experience, levels=c("Has previously served in the government","Has not previously served in the government")))

# Subset to just always responders
tun_bjka_econ_long_always <- tun_bjka_econ_long_na2[tun_bjka_econ_long_na2$cjoint_response=="Always",]


library(cregg)

# save formula
f1 <- choice ~ Sex + Policy + Age + Education + Job + Experience + Town + Party

# mm for patriarchal vs. egalitarian
mm_patr_econ <- cj(tun_bjka_econ_long_always, f1, id = ~ responseid, estimate = "mm", by = ~ type, level_order="descending")

mm_diff_patr_econ <- cj(tun_bjka_econ_long_always, f1, id = ~responseid, estimate = "mm_diff", 
                   by = ~type, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_1 <- mm_patr_econ[plot_vars]
mm_2 <- mm_diff_patr_econ[plot_vars]
mm_econ_all <- rbind(mm_1,mm_2)

mm_econ_all$BY <- factor(mm_econ_all$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_econ_all$BY), vl=c(0.5,0.5,0.0)) 

my_breaks <- function(q) { if (max(q) < .35) seq(-.1, .2, .1) else seq(.3, .7, .1) }

mm_plot_econ<-plot(mm_econ_all,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                           size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  #scale_x_continuous(breaks = scales::pretty_breaks(5))
  scale_x_continuous(breaks = my_breaks) +
  ggplot2::theme(axis.text.x=element_text(size=7),
                 axis.title.x=element_text(size=11),
                 axis.text.y=element_text(size=8),
                 axis.title.y=element_text(size=11),
                 panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 panel.background = element_blank(),
                 legend.position="bottom",
                 legend.title = element_blank())

mm_plot_econ

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_mm_econ.png",
    width = 8, height = 6, units = 'in', res = 300)
mm_plot_econ
dev.off()


#########################################################################   
#     ROBUSTNESS CHECK 2: USE RESIDUALS FOR ALTERNATIVE MEASURE         #   
#########################################################################

#########################################################################   
#             remove NAs for robustness check (BJKA only)               #   
#########################################################################

tun_bjka <- tun_bjka %>% 
  filter(!is.na(resp_fem)) %>% 
  filter(!is.na(age)) %>%
  filter(!is.na(educ_high)) %>% 
  filter(!is.na(resp_mar_wid)) %>%
  filter(!is.na(vote_yes)) %>% 
  filter(!is.na(employ_status)) %>%
  filter(!is.na(patr_index))


#########################################################################   
#                 use residuals to define types                         #   
#########################################################################

# Controlling for education, marital status, voting, employment
patr_bjka2 <- lm(patr_index ~ resp_fem + age + educ_high +
                   resp_mar_wid + vote_yes + employ_status + del_en, data = tun_bjka)
summary(patr_bjka2)

# Controlling for education, marital status, voting, employment
patr_yg2 <- lm(patr_index ~ resp_fem + age + educ_high +
                 resp_mar_wid + vote_yes + employ_status + city_en, data = tun_yg)
summary(patr_yg2)


tun_bjka$predicted <- predict(patr_bjka2)   # Save the predicted values
tun_bjka$residuals <- residuals(patr_bjka2) # Save the residual values

tun_yg$predicted <- predict(patr_yg2)   # Save the predicted values
tun_yg$residuals <- residuals(patr_yg2) # Save the residual values

# Quick look at the actual, predicted, and residual values
tun_bjka %>% dplyr::select(patr_index, predicted, residuals) %>% head()
tun_yg %>% dplyr::select(patr_index, predicted, residuals) %>% head()

median(tun_bjka$residuals)
median(tun_yg$residuals)

tun_bjka$patr_resid <- ifelse(tun_bjka$residuals >= median(tun_bjka$residuals),1,0)
tun_bjka$type_resid <- ifelse(tun_bjka$patr_resid==1, "Patriarchal","Egalitarian")

tun_yg$patr_resid <- ifelse(tun_yg$residuals >= median(tun_yg$residuals),1,0)
tun_yg$type_resid <- ifelse(tun_yg$patr_resid==1, "Patriarchal","Egalitarian")


g2 <- ggplot(tun_bjka) +
  geom_hline(yintercept = 0,
             color = "gray70",
             linetype = "dashed") +
  geom_point(aes(x = patr_index,
                 y = residuals,
                 color = factor(patr_resid, 
                                label = c("Egalitarian", "Patriarchal")),
                 shape = factor(patr_resid, 
                                label = c("Egalitarian", "Patriarchal"))),
             alpha = 0.8,
             na.rm = TRUE) +
  scale_color_manual(values = c("black", "darkgray")) +
  theme_few() +
  theme(legend.position = c(0.8, 0.2)) +
  labs(x = "Gender Egalitarian Index",
       y = "OLS residuals",
       color = "Type", 
       shape = "Type")

g2

g3 <- ggplot(tun_yg) +
  geom_hline(yintercept = 0,
             color = "gray70",
             linetype = "dashed") +
  geom_point(aes(x = patr_index,
                 y = residuals,
                 color = factor(patr_resid, 
                                label = c("Egalitarian", "Patriarchal")),
                 shape = factor(patr_resid, 
                                label = c("Egalitarian", "Patriarchal"))),
             alpha = 0.8,
             na.rm = TRUE) +
  scale_color_manual(values = c("black", "darkgray")) +
  theme_few() +
  theme(legend.position = c(0.8, 0.2)) +
  labs(x = "Gender Egalitarian Index",
       y = "OLS residuals",
       color = "Type", 
       shape = "Type")

g3


#########################################################################   
#       robustess check 2: BJKA vignette analysis                       #   
#########################################################################

m_all <- lm(support ~ tassign_women, data=tun_bjka)
summary(m_all)

m_patr <- lm(support ~ tassign_women, data=tun_bjka[tun_bjka$patr_resid==1,])
summary(m_patr)

m_egal <- lm(support ~ tassign_women, data=tun_bjka[tun_bjka$patr_resid==0,])
summary(m_egal)


#########################################################################   
#       robustess check 2: BJKA conjoint analysis                       #   
#########################################################################

tun_bjka$never_resp <- ifelse(is.na(tun_bjka$choice_1)==T &
                                   is.na(tun_bjka$choice_3)==T &
                                   is.na(tun_bjka$choice_5)==T &
                                   is.na(tun_bjka$choice_7)==T, 1,0)


tun_bjka$sometimes_resp <- ifelse(is.na(tun_bjka$choice_1)==T |
                                       is.na(tun_bjka$choice_3)==T |
                                       is.na(tun_bjka$choice_5)==T |
                                       is.na(tun_bjka$choice_7)==T, 1,0)

tun_bjka$cjoint_response <- ifelse(tun_bjka$sometimes_resp==1,"Sometimes","Always")
tun_bjka$cjoint_response <- ifelse(tun_bjka$never_resp==1,"Never",tun_bjka$cjoint_response)
table(tun_bjka$cjoint_response)


tun_bjka_long <- reshape(data = tun_bjka, 
                            idvar = "responseid", 
                            varying = list(Sex=c("sex_1","sex_2","sex_3","sex_4","sex_5","sex_6","sex_7","sex_8"),
                                           Party=c("party_1","party_2","party_3","party_4","party_5","party_6","party_7","party_8"),
                                           Job=c("job_1","job_2","job_3","job_4","job_5","job_6","job_7","job_8"),
                                           Town=c("town_1","town_2","town_3","town_4","town_5","town_6","town_7","town_8"),
                                           Education=c("educ_1","educ_2","educ_3","educ_4","educ_5","educ_6","educ_7","educ_8"),
                                           Age=c("age_1","age_2","age_3","age_4","age_5","age_6","age_7","age_8"),
                                           Policy=c("policy_1","policy_2","policy_3","policy_4","policy_5","policy_6","policy_7","policy_8"),
                                           Experience=c("exper_1","exper_2","exper_3","exper_4","exper_5","exper_6","exper_7","exper_8"),
                                           choice=c("choice_1","choice_2","choice_3","choice_4","choice_5","choice_6","choice_7","choice_8"),
                                           contest=c("contest_1","contest_2","contest_3","contest_4","contest_5","contest_6","contest_7","contest_8")), 
                            direction="long", 
                            v.names = c("Sex","Party","Job","Town","Education","Age","Policy","Experience","choice","contest"),
                            sep="_")

tun_bjka_long_na2 <- tun_bjka_long %>% filter(!is.na(choice)) %>% filter(!is.na(patriarchal)) %>%
  mutate(Sex = factor(Sex, levels=c("Male","Female"))) %>%
  mutate(Party = factor(Party, levels=c("Initiative Party","Afek Tounes","UPL","Democratic Current",
                                        "Current of Love","al-Irada Movement","Popular Front",
                                        "Machrouu Tounes","Ennahda","Nidaa Tounes"))) %>%
  mutate(Job = factor(Job, levels = c("Director", "Lawyer","Teacher","Employee in the private sector","Political Activist","Union Activist"))) %>%
  mutate(Town = factor(Town, levels = c("within 10 km of here","within 20 km of here","within 30 km of here","within 40 km of here"))) %>%
  mutate(Education = factor(Education, levels = c("Bachelors","Masters","Studied outside of Tunisia"))) %>%
  mutate(Age = factor(Age, levels=c("70","60","50","40","30"))) %>%
  mutate(Policy = factor(Policy, levels=c("Improve security","Improve women's rights","Fight corruption","Improve the economic situation"))) %>%
  mutate(Experience = factor(Experience, levels=c("Has previously served in the government","Has not previously served in the government")))

# Subset to just always responders
tun_bjka_long_always <- tun_bjka_long_na2[tun_bjka_long_na2$cjoint_response=="Always",]


library(cregg)

# save formula
f1 <- choice ~ Sex + Policy + Age + Education + Job + Experience + Town + Party

# mm for patriarchal vs. egalitarian
mm_patr <- cj(tun_bjka_long_always, f1, id = ~ responseid, estimate = "mm", by = ~ type_resid, level_order="descending")
plot(mm_patr, vline = 0.5)

mm_diff_patr <- cj(tun_bjka_long_always, f1, id = ~responseid, estimate = "mm_diff", 
                   by = ~type_resid, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_1 <- mm_patr[plot_vars]
mm_2 <- mm_diff_patr[plot_vars]
mm_bjka_all <- rbind(mm_1,mm_2)

mm_bjka_all$BY <- factor(mm_bjka_all$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_bjka_all$BY), vl=c(0.5,0.5,0.0)) 

my_breaks <- function(q) { if (max(q) < .3) seq(-.1, .2, .1) else seq(.3, .7, .1) }

mm_plot_bjka<-plot(mm_bjka_all,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                           size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  #scale_x_continuous(breaks = scales::pretty_breaks(5))
  scale_x_continuous(breaks = my_breaks)+
  ggplot2::theme(axis.text.x=element_text(size=7),
                 axis.title.x=element_text(size=11),
                 axis.text.y=element_text(size=8),
                 axis.title.y=element_text(size=11),
                 panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 panel.background = element_blank(),
                 legend.position="bottom",
                 legend.title = element_blank())

mm_plot_bjka

png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_mm_resid_bjka.png",
    width = 8, height = 6, units = 'in', res = 300)
mm_plot_bjka
dev.off()


#########################################################################   
#       robustess check 2: YouGov vignette analysis                     #   
#########################################################################

tun_yg_long <- gather(tun_yg, condition, support, woman_control:man_security, factor_key=TRUE)

m_all <- lm(support ~ condition + respondentid, data=tun_yg_long)
summary(m_all)

m_patr <- lm(support ~ condition + respondentid, data=tun_yg_long[tun_yg_long$patr_resid==1,])
summary(m_patr)

m_egal <- lm(support ~ condition + respondentid, data=tun_yg_long[tun_yg_long$patr_resid==0,])
summary(m_egal)

# Cluster by respondent
overall_clust <- super.cluster.fun(m_all, tun_yg_long$respondentid)
patr_clust <- super.cluster.fun(m_patr, tun_yg_long$respondentid[tun_yg_long$patr_resid==1])
egal_clust <- super.cluster.fun(m_egal, tun_yg_long$respondentid[tun_yg_long$patr_resid==0])

stargazer(overall_clust,patr_clust,egal_clust,
          title = "Vignette (YouGov, 2019)",
          dep.var.labels   = c("Support for Candidate"),
          omit = c("respondentid"),
          notes = "SEs clustered by respondent",
          add.lines = list(c("Respondent FE" ,"Yes","Yes", "Yes", "Yes")))

#########################################################################   
#       robustess check 2: YouGov conjoint analysis                     #   
#########################################################################

tun_yg$choice_profileA <- ifelse(tun_yg$QS3_cj1=="Candidate A",1,0)
tun_yg$choice_profileB <- ifelse(tun_yg$QS3_cj1=="Candidate B",1,0)

tun_yg$choice_profileC <- ifelse(tun_yg$QS3_cj2=="Candidate C",1,0)
tun_yg$choice_profileD <- ifelse(tun_yg$QS3_cj2=="Candidate D",1,0)

tun_yg$choice_profileE <- ifelse(tun_yg$QS3_cj3=="Candidate E",1,0)
tun_yg$choice_profileF <- ifelse(tun_yg$QS3_cj3=="Candidate F",1,0)

tun_yg$choice_profileG <- ifelse(tun_yg$QS3_cj4=="Candidate G",1,0)
tun_yg$choice_profileH <- ifelse(tun_yg$QS3_cj4=="Candidate H",1,0)

tun_yg$choice_profileI <- ifelse(tun_yg$QS3_cj5=="Candidate I",1,0)
tun_yg$choice_profileJ <- ifelse(tun_yg$QS3_cj5=="Candidate J",1,0)


tun_yg_conjoint <- reshape(data = tun_yg,
                           idvar = "RecordNo", 
                           varying = list(Sex=c("A1_profileA","A1_profileB","A1_profileC","A1_profileD","A1_profileE","A1_profileF",
                                                "A1_profileG","A1_profileH","A1_profileI","A1_profileJ"),
                                          Party=c("A2_profileA","A2_profileB","A2_profileC","A2_profileD","A2_profileE","A2_profileF",
                                                  "A2_profileG","A2_profileH","A2_profileI","A2_profileJ"),
                                          Job=c("A3_profileA","A3_profileB","A3_profileC","A3_profileD","A3_profileE","A3_profileF",
                                                "A3_profileG","A3_profileH","A3_profileI","A3_profileJ"),
                                          Town=c("A4_profileA","A4_profileB","A4_profileC","A4_profileD","A4_profileE","A4_profileF",
                                                 "A4_profileG","A4_profileH","A4_profileI","A4_profileJ"),
                                          Education=c("A5_profileA","A5_profileB","A5_profileC","A5_profileD","A5_profileE","A5_profileF",
                                                      "A5_profileG","A5_profileH","A5_profileI","A5_profileJ"),
                                          Experience=c("A6_profileA","A6_profileB","A6_profileC","A6_profileD","A6_profileE","A6_profileF",
                                                       "A6_profileG","A6_profileH","A6_profileI","A6_profileJ"),
                                          Policy=c("A7_profileA","A7_profileB","A7_profileC","A7_profileD","A7_profileE","A7_profileF",
                                                   "A7_profileG","A7_profileH","A7_profileI","A7_profileJ"),
                                          Age=c("A8_profileA","A8_profileB","A8_profileC","A8_profileD","A8_profileE","A8_profileF",
                                                "A8_profileG","A8_profileH","A8_profileI","A8_profileJ"),
                                          choice=c("choice_profileA","choice_profileB","choice_profileC","choice_profileD","choice_profileE",
                                                   "choice_profileF","choice_profileG","choice_profileH","choice_profileI","choice_profileJ")), 
                           direction="long", 
                           v.names = c("Sex","Party","Job","Town","Education","Experience","Policy","Age","choice"),
                           sep="_")

# Rename Jabha Shabiyya --> Popular Front
tun_yg_conjoint$Party <- ifelse(tun_yg_conjoint$Party=="Jabha Shabiyya", "Popular Front",
                                tun_yg_conjoint$Party)

# Rename Studied outside of Tunis --> Studied outside of Tunisia
tun_yg_conjoint$Education <- ifelse(tun_yg_conjoint$Education=="Studied outside of Tunis", "Studied outside of Tunisia",
                                    tun_yg_conjoint$Education)

# Rename Director of an organization --> Director
tun_yg_conjoint$Job <- ifelse(tun_yg_conjoint$Job=="Director of an organization", "Director",
                              tun_yg_conjoint$Job)

tun_yg_conjoint_na <- tun_yg_conjoint %>% filter(!is.na(choice)) %>% filter(!is.na(patriarchal)) %>%
  mutate(Sex = factor(Sex, levels=c("Male","Female"))) %>%
  mutate(Party = factor(Party, levels=c("Independent List","Democratic Current","Popular Front",
                                        "Ennahda","Nidaa Tounes"))) %>%
  mutate(Job = factor(Job, levels = c("Director", "Lawyer","Teacher",
                                      "Private sector employee","Political Activist","Union Activist"))) %>%
  mutate(Town = factor(Town, levels = c("Within 10 km of your municipality",
                                        "Within 20 km of your municipality",
                                        "Within 30 km of your municipality",
                                        "Within 40 km of your municipality"))) %>%
  mutate(Education = factor(Education, levels = c("Bachelors","Masters","Studied outside of Tunisia"))) %>%
  mutate(Age = factor(Age, levels=c("71","60","51","42","33"))) %>%
  mutate(Policy = factor(Policy, levels=c("Improve security","Improve women's rights","Fight corruption","Improve the economic situation"))) %>%
  mutate(Experience = factor(Experience, levels=c("Has previously served in the government","Has not previously served in the government")))


library(cregg)

# save formula
f1 <- choice ~ Sex + Policy + Age + Education + Job + Experience + Town + Party

# mm for patriarchal vs. egalitarian
mm_patr <- cj(tun_yg_conjoint_na, f1, id = ~ RecordNo, estimate = "mm", by = ~ type_resid, level_order="descending")
plot(mm_patr, vline = 0.5)

mm_diff_patr <- cj(tun_yg_conjoint_na, f1, id = ~RecordNo, estimate = "mm_diff", 
                   by = ~type_resid, level_order="descending")

plot_vars <- c("BY","statistic","feature", "level","estimate", "std.error",
               "z","p","lower","upper")

mm_1 <- mm_patr[plot_vars]
mm_2 <- mm_diff_patr[plot_vars]
mm_yg_all <- rbind(mm_1,mm_2)

mm_yg_all$BY <- factor(mm_yg_all$BY,levels=c("Patriarchal","Egalitarian","Patriarchal - Egalitarian"))
vline.dat <- data.frame(BY=levels(mm_yg_all$BY), vl=c(0.5,0.5,0.0)) 

my_breaks <- function(q) { if (max(q) < .3) seq(-.1, .2, .1) else seq(.3, .7, .1) }

mm_plot_yg<-plot(mm_yg_all,vline=NULL) + ggplot2::geom_vline(aes(xintercept=vl), data=vline.dat,
                                                           size=.25, linetype='dashed') + 
  ggplot2::facet_wrap(~BY, ncol = 3L, scales = "free_x")  +
  #scale_x_continuous(breaks = scales::pretty_breaks(5))
  scale_x_continuous(breaks = my_breaks)+
  ggplot2::theme(axis.text.x=element_text(size=7),
                 axis.title.x=element_text(size=11),
                 axis.text.y=element_text(size=8),
                 axis.title.y=element_text(size=11),
                 panel.grid.major = element_blank(),
                 panel.grid.minor = element_blank(),
                 panel.background = element_blank(),
                 legend.position="bottom",
                 legend.title = element_blank())

mm_plot_yg


png(file = "~/Dropbox/Voters and Female Politicians/Political Behavior/Figures/conjoint_mm_resid_yg.png",
    width = 8, height = 6, units = 'in', res = 300)
mm_plot_yg
dev.off()






